Screening and identification of potential key biomarkers for glucocorticoid-induced osteonecrosis of the femoral head

Background Glucocorticoid-induced osteonecrosis of the femoral head (GIONFH) is a common disease in osteoarticular surgery, with a high disability rate, which brings great physical and mental pain and economic burden to patients. Its specific pathogenesis has not been fully demonstrated, and there is a lack of recognized effective biomarkers for earlier detection and prompt treatment. This has become an urgent clinical problem for orthopedic scholars. Materials and methods We downloaded the gene expression profile dataset GSE123568 from the Gene Expression Omnibus database, used STRING and Cytoscape to carry out module analysis and built a gene interaction network. The four core genes most related to GIONFH in this network were ultimately found out by precise analysis and animal experiment were then conducted for verification. In this verification process, thirty-six New Zealand white rabbits were randomly divided into blank control group, model group and drug group. Except for the blank control group, the animal model of GIONFH was established by lipopolysaccharide and methylprednisolone, while the drug group was given the lipid-lowering drugs for intervention as planned. The rabbits were taken for magnetic resonance imaging at different stages, and their femoral head specimens were taken for pathological examination, then the expression of target genes in the femoral head specimens of corresponding groups was detected. Validation methods included RT-PCR and pathological examination. Results A total of 679 differential genes were selected at first, including 276 up-regulated genes and 403 down-regulated genes. Finally, four genes with the highest degree of correlation were screened. Animal experiment results showed that ASXL1 and BNIP3L were in low expression, while FCGR2A and TYROBP were highly expressed. Conclusion Through animal experiments, it was confirmed that ASXL1, BNIP3L, FCGR2A and TYROBP screened from the comparative analysis of multiple genes in the database were closely related to GIONFH, which is important for early diagnosis of Glucocorticoid-induced osteonecrosis of the femoral head.


Introduction
Glucocorticoids, a form of steroid hormones, are mainly used in the treatment of autoimmune diseases, such as rheumatoid arthritis, systemic lupus erythematosus and scleroderma, as well as organ transplantation and severe infection. Long-term or extensive use of glucocorticoids will bring more complications to patients, among which GIONFH is an example [1][2][3][4]. As a common type of nontraumatic osteonecrosis, GIONFH may result in the collapse of the subchondral bone structure without immediate treatment. Previous studies have shown that steroid hormones therapy can increase the risk of osteonecrosis by 4-50% [5].
In recent years, more attention has been gradually paid to GIONFH [6][7][8]. In the fight against Severe Acute Respiratory Syndrome (SARS) in 2003, glucocorticoids were used at high doses to save a large number of critically ill patients. Unfortunately, these steroid hormones also had side effects-the incidence of GIONFH was as high as 24-30% [9][10][11]. As for the novel coronavirus disease 2019 (COVID-19) epidemic, which is similar to SARS, steroid hormone is still an important medication for the treatment of infections to rescue lives, who will also be at risk of developing GIONFH in the future [12][13][14]. Therefore, early diagnosis and timely prevention and treatment of GIONFH have important clinical significance.
However, the pathogenesis of GIONFH is yet to be seen. The hypotheses proposed by researchers mainly include as follows-intravascular coagulation, vascular endothelial injury, increased intraosseous pressure, lipid metabolism disorder etc. [15]. Among them, the disorder of lipid metabolism is an important manifestation in GIONFH [16,17]. It was confirmed in previous studies that steroid use may result in changes in the lipid metabolism, which increase the risk of fat embolism [18,19]. Our team's studies also confirmed that GIONFH is closely related to the pathogenesis of lipid metabolism disorder induced by high-dose hormone in vivo [20][21][22]. Therefore, in order to further confirm the association between the onset of GIONFH and lipid metabolism disorder, a hypolipidemic drug group was added as comparative research to the Glucocorticoid-induced femoral head necrosis model rabbits in the animal experiment of this study.
In the current clinical research, there is a lack of methods and standards for early diagnosis of GIONFH. It is asymptomatic during the early phases. The first symptom is pain in the hip, which tends to be in proportion with lesion size, and in general precedes the beginning of femoral head collapse, which occurs on average after 8 months [23,24]. At the moment of hip pain and dysfunction, and obvious osteonecrosis and collapse of the femoral head appear in pathological and imaging examinations, in which case surgical treatment is badly needed [25]. With the wide application of bioinformatics and the continuous development of microarray technology at the genome level, researchers have sequenced more and more disease samples with increasing fragment information, and the public functional genomics database has been greatly enriched [26,27]. By screening and mining the massive information, we can find the mechanism of the occurrence and development of non-tumor diseases from the gene level, which helps us to screen out potential biomarkers. Then we verified the key markers screened before through animal experiments, it is expected to find new methods to facilitate the early diagnosis, detection and prevention of GIONFH.

Materials and methods
Gene chip data GEO (http:// www. ncbi. nlm. nih. gov/ geo) is a public functional genomics database with high storage throughout gene expression data, chips, and microarrays [28]. The gene expression profile dataset GSE123568 was downloaded from the GEO database, and probes were converted to corresponding gene symbols based on annotation information on the platform [29]. The GSE123568 dataset contained tissue specimens from 30 GIONFH patients and 10 non-GIONFH patients. GEO2R (http:// www. ncbi. nlm. nih. gov/ geo/ geo2r), an interactive web tool, was used to screen DEGs between GIONFH and non-GIONFH samples. GEO2R allows users to compare two or more datasets in a GEO series in order to identify DEGs under certain conditions. Probe sets without corresponding gene symbols or genes with more than one probe set were removed, respectively. A log FC (fold change) > 1 or log FC < − 1 with an adjusted P value < 0.01 was considered statistically significant.

PPI network and module analysis
STRING is a database for predicting known protein interactions. Its interactions include direct (physical) associations and indirect (functional) associations. Analysis of function and interactions between proteins with STRING may provide insights into the mechanisms involved in disease development or progression [30]. In this study, the PPI network of DEGs was constructed using the STRING database, and interactions with a composite score (the support of the data) > 0.4 were considered statistically significant. Cytoscape (3.8.0 version) is an open-source software platform for visualizing complex networks and integrating them with any type of attribute data, including bioinformatics and social network analysis [31]. Cytoscape's plugin MCODE (2.0.0 version), is used to cluster a given network based on topology to find densely connected regions (highly interconnected regions) [32]. The PPI network as well as the biological process map were drawn with cytoscape, and the most important modules in the PPI network were identified with MCODE. The selection criteria were as follows: MCODE score > 5, degree cutoff = 2, node score cutoff = 0.2, maximum depth = 100, and k score = 2. Biological process analysis of central genes and protein network interaction mapping were performed by BiNGO, a plugin of cytoscape. The relationships among the nine hub genes were explored by STRING.

GO and KEGG enrichment analysis of DEGs
DAVIDv6.8 (http: //david.ncifcrf.gov/) [33], a visual integrated discovery database with annotation contains a complete update of knowledge base. KEGG is a database used for understanding advanced functions and practical experimental techniques in biological systems (e.g., cells, organisms and ecosystems) from molecular-level information (especially large-scale molecular datasets generated by genome sequencing and other high-throughput processes) [34]. The Gene Ontology (GO) knowledge base has the maximum information about gene function around the world, which can be recognized both by human and machine. It is also the basis for computational analysis of large-scale molecular biology and genetic experiments in biomedical research [35]. Functional enrichment analysis of DEGs was performed using the DAVID online database, and P < 0.05 was considered statistically significant.

Selection and analysis of key genes
UniProt (universal protein) is a database (https:// www. unipr ot. org/) collecting protein resources and interconnected with other databases. It has one of the most extensive protein sequences with most comprehensive annotations so far. The 18 genes we screened were consulted using the UniProt database and annotated in lists. Functional enrichment analysis of potential genes was then performed by DAVID [28] based on up-regulated genes (EPB41, BNIP3L, KLF1, SLC7A5, C10orf10, ASXL1) and down-regulated genes (CD86, FCGR3A, FCGR2B, CYBB, TLR2, TLR4, ITGAX, TLR8, TLR1, FCGR2A, TYROBP, MS4A6A). Such analysis was also performed for individual genes (BNIP3L, ASXL1, FCGR2A, TYROBP) to search for possible protein pathways.

Animal experiments
Thirty-six adult New Zealand white rabbits weighing 2.5-3.0 kg were purchased from the Animal Experimental Center of Jiangxi University of Traditional Chinese Medicine with license number of SYXK[gan]2004-0001. The 18 male and 18 female rabbits were single-housed with constant temperature at 20-23℃ under the humidity of 60%, and the light-dark cycle was 12 h each. After one week of adaptive feeding, those rabbits were randomly divided into three groups: blank control group (n = 12), model group (n = 12) and drug group (n = 12), half male and half female in each groups. Rabbits in the model group and drug group were intravenously injected with lipopolysaccharide (LPS, Sigma, USA) (10 μg/ kg) via the marginal ear vein, and the second dose was administered 24 h later. After two injections of LPS, methylprednisolone sodium succinate (MSS, Huazhong Pharmaceutical, China) was injected intramuscularly at a dose of 20 mg/kg with a 24 h interval between each injection for a total of three injections. During hormone injection, 80,000 U of penicillin was intramuscularly injected into the right buttock to prevent infection (3 days) [36,37]. After the completion of hormone injection, the drug group was intragastrical administered with simvastatin (Jingxin Pharmaceutical, China) water solution every day at a dose of 2 mg/10 ml per time for 4 weeks, the control and model groups were intragastrical administered with the same amount of saline. Magnetic resonance images of the hip joints and femoral head specimens were collected from three groups of rabbits at different time (after 2, 4, and 8 weeks). This study was conducted according to the recommendations of the National Institutes of Health Guide for the Care and Use of Laboratory Animals. The protocol was approved by the Experimental Animal Ethics Committee of the Fourth Affiliated Hospital of Nanchang University.

Magnetic resonance imaging scanning
Magnetic resonance imaging (MRI) examination was performed at two, four, and eight weeks after the first hormone injection. 3% pentobarbital sodium solution (35 mg/kg) was injected intravenously into the ear margin of rabbits. After the anesthesia took effect, the bilateral femoral heads of all rabbits were scanned by a 3.0 T magnetic resonance imaging (MRI) scanner (Achieva TX, Philips). The parameters were as follows: T1WI (TR 700 ms, TE 45 ms), T2WI (TR 2600 ms, TE 200 ms), Stir (TR 1500 ms, TE 60 ms), slice thickness 1.2 mm, flip angle 90° and matrix (reconstruction 256 × matrix scan 256 mm). Then the MRI results of GIONFH in each group were compared.

Extraction and anatomy of specimens
The rabbits were sacrificed by air embolism at two, four, and eight weeks after hormone injection, and the femoral heads were harvested. Washed with ice phosphate buffer saline (PBS), the left femoral heads of all rabbits were preserved in liquid nitrogen for quantitative real-time PCR. Another femoral head was collected and immediately

Pathological examination
After the femoral heads were fixed with 4% paraformaldehyde for 24 h, they were washed with 0.2 M PBS (pH 7.4) again and then decalcified in 10% ethylene diamine tetraacetic acid (EDTA) and neutralized with sodium sulphate buffer for approximately 30 days. After decalcification, the tissues were embedded in paraffin, cut into thick sections of 5 mm, and stained with hematoxylin and eosin. The osteonecrosis of the femoral head in H&E staining specimens was observed under a microscope, and the microphotos were taken. Leica Co. W550CW signal acquisition and analysis system (Weztlar, Germany) was used for observation, analysis and evaluation.

Real-time PCR
Total RNA was extracted using rnaiso plus (Takara bio, Japan) according to the manufacturer's instructions. The RNA concentration was evaluated by A260/A280 measurement, then reverse transcription was performed, and the reverse transcription product was amplified by qPCR Expression of BNIP3L, ASXL1, FCGR2A and TYROBP were analyzed using SYBR premix Ex TaqTM II (Takara bio, Japan) on ABI 7500 rapid real-time PCR system (applied biological systems, USA). Obtain the cycle threshold (CT) value and use 2-ΔΔ the relative expression of each target gene was calculated by CT method, and the data were normalized to GAPDH expression. The primer sequence is shown in Table 1.

Statistical analysis
Quantitative data are presented as mean ± standard deviation (SD). Hardy-Weinberg genetic balance method to test population representativeness, the phenotypic differences between groups were analyzed by paired t-test to estimate the relationship between each genotype and the risk of GIONFH the odds ratio (OR), 95% confidence interval (CI) and P value were calculated. Data were analyzed using SPSS 23.0 software (SPSS, USA). A level of P < 0.05 was considered to be significant.

Analysis of gene arrays
Through GEO2R analysis, 679 DEGs were selected including 276 up-regulated genes and 403 down-regulated genes; the volcano plot (Fig. 1A) shows the relationship between fold change and P value in this set of microarray groups. The 10 genes with the highest weights were all down-regulated genes by Cytoscape analysis, which were CD86, FCGR3A, FCGR2B, CYBB, TLR2, TLR4, ITGAX, TLR8, FCGR2A, and TLR1 (Fig. 1B); 9 groups were found by MCODE, and there was a hub gene in each group, of which the up-regulated genes were ASXL1, BNIP3L, C10orf10, EPB41, KLF1, and SLC7A5, and the down-regulated genes were FCGR2A, MS4A6A, and TYROBP ( Fig. 2A).

Functional enrichment analysis
It should be noted that the above two groups of genes have an overlapping gene (FCGR2A), so a total of 18 differential genes are synthesized from the two groups. When exploring the relationship between the nine hub genes, some genes were found to be interrelated (Fig. 2B). A protein interaction network was drawn for the 18 differential genes to provide protein expression or potential pathways and targets that may be associated with GIONFH (Fig. 3). The results of functional enrichment analysis of 18 genes mainly focus on down-regulated genes (Table 1), of which important parts have been marked (Fig. 4). The results showed that eighteen differential genes were annotated, among which BNIP3L, ASXL1, FCGR2A and TYROBP were highly correlated with osteonecrosis (Tables 2 and 3).

Magnetic resonance imaging results
Magnetic resonance imaging (MRI) is considered the most sensitive, specific modality and golden standard for the diagnosis and evaluation of GIONFH. It was used to analyze the imaging changes in the rabbits' femoral heads. The results of T2 weighted MRI showed that in the blank control group, the shape of bilateral femoral heads was regular, the surface was smooth, no abnormal signals were found, and the bilateral hip joint space was normal ( Fig. 5 A, B). As for the model group after 8 weeks, the T2-weighted images showed spotty, heterogeneous high signal intensity, slight collapse and flattening of bilateral femoral heads, edema signals in bilateral joint cavities (red arrow), and uneven joint spaces (Fig. 5 C, D). In the drug group, the joint space were normal, no obvious effusion was in the hip cavity, the epiphyseal line of the femoral heads were slightly blurred, and no obvious abnormal signals were found (Fig. 5 E, F). These results indicated that steroid hormone-induced osteonecrosis of the femoral head, however, simvastatin could exert a protective effect on GIONFH in vivo.

Histopathological analysis H&E staining
Based on the diagnostic criteria of GIONFH, we evaluated femoral head specimens of rabbits using H&E staining. In this study, none of the rabbits in the control group developed osteonecrosis. The bone trabeculas were regularly arrayed, with complete structure, clearly visible osteocytes and a few empty lacunaes (Fig. 6 A, B).
In the model group, typical osteonecrosis occurred. The lipocytes were enlarged. Bone trabeculas turned thinner, and many empty lacunaes were observed (Fig. 6 C,  D). In the drug group, there were a few lipocytes, and no apparent necrotic debris was noted. The bone trabeculas regularly arrayed, and only a few empty lacunaes were observed (Fig. 6 E, F).  Expression of the key differential genes in the rabbit model of GIONFH In order to investigate the expression of the key differential genes related to GIONFH in rabbits, real-time PCR was performed to measure the expression level of the ASXL1, BNIP3L, FCGR2A and TYROBP at different time points (after 2, 4, and 8 weeks). After 2 weeks, the results showed that the expression level of ASXL1 and BNIP3L were significantly decreased and the TYROBP robust increased in the model group compared to the blank control group (Fig. 7 A). But this effect was abolished when the rabbits were treated with simvastatin in the drug group. Similarly, after 4 and 8 weeks, the expression level of ASXL1 and BNIP3L decreased significantly in the model group compared to the blank control group, while the FCGR2A and TYROBP increased remarkable (Fig. 7  B, C). It should be noted that this effect was blocked by simvastatin in the drug group. Together, these data suggested that the four key differential genes ASXL1, BNIP3L, FCGR2A and TYROBP were significantly associated with osteonecrosis of the femoral head.

Discussion
Although steroid hormones play an important role in the treatment of many diseases, they produce a series of side effects. For example, long-term or extensive use of steroid hormone may cause GIONFH [38]. Previous studies have shown that without early intervention, 80% of patients will develop a collapsed femoral head and significantly impaired hip function within 1-5 years of onset [39]. Such treatments for collapsed and necrotic femoral head as core decompression, Osteotomies, bone pedicled vascularized or autologous stem cell implantation are all invasive operations [40][41][42][43]. After the operations patients are very likely to have unbalanced bone resorption, which, as is well known, can be enhanced by corticosteroids [42]. Due to the negative effect in the long run, total hip arthroplasty (THA) is eventually an inevitable option. Therefore, more and more researchers are devoted to the study of early diagnosis and preventive drugs for GIONFH. Bioinformatics analysis and microarray technology at the genome level are mostly adopted in tumor research. At first, the researchers sequenced the tumor sample cells of human, and uploaded the sequenced gene array to the gene bank. Then, a variety of sequencing results were downloaded from the gene bank by other researchers and tens of thousands of gene fragments were obtained after researchers' classification and annotation. At last, key genes or mutated genes were identified from these gene fragments for research to discover the disease pathogenic mechanism and therapeutic targets [44,45]. With the help of the analysis platform of lipid metabolomics technology, our research group has studied the variation trend of the metabolic profile spectrum of animal samples of GIONFH before and after the intervention, and confirmed that GIONFH is closely related to the expression of lipid metabolism gene. Moreover, animal experimental studies have confirmed that the application of lipid-lowering drugs could effectively improve the occurrence and disease progression of GIONFH [20,21]. In this study, we obtained 679 differentially expressed genes by screening between the gene array of human GIONFH and non-GIONFH. In the screening analysis of 679 differentially expressed genes, the following two algorithm were used to identify the potential key genes. Ten weight genes were obtained by MCC (maximum clique centrality) algorithm. According to the statistics on its official website, the key networks captured by MCC algorithm will bring new insights into the basic regulatory networks and protein drug targets for experimental biologists. In addition, we also used the MCODE (molecular complex detection) Table 3, we selected some important biological processes, cell structure, molecular functions and signal pathways from the functional enrichment analysis. A Pie chart is the relationship between items and FDR. B The relationship between bar graph entries and P-values  plug-in to obtain nine gene groups. Despite the complexity of running the code of MCODE, this method has been used by many medical researchers for many times, which is enough to verify its reliability. Therefore, we used this method to obtain 18 potential genes to be verified. We compared and analyzed the correlation between the related proteins expressed by the 18 genes and steroid-induced femoral head necrosis one by one. TYROBP and FCGR2A were associated with osteonecrosis at different sites in the same pathway, among which the mutation of TYROBP-DAP12 encodes membrane receptor component cells in natural killer cells and myeloid cells, which has been confirmed by Juha Paloneva et al. [46], meaning TYROBP-mediated signaling pathway played an vital role in human bone tissue. Another study by Juha Paloneva and her team [47] also demonstrated an important role for the DAP12-TREM2 signaling complex in the differentiation and function of osteoclasts. This undoubtedly highlights the importance of the high expression of TYROBP in the rabbit GIONFH model tissue. FCGR2A acts as an Fc region-binding receptor of immunoglobulin γ and has been found to be associated with many inflammatory diseases [48]. FcgammaRIIa (FcγRII) is a multivalent IgG receptor expressed primarily by myeloid cells, and its binding to lipid rafts (microdomains rich in cholesterol and sphingolipids) is critical for efficient signaling of the pathway [49]. Yang et al. [50] found in the study of myeloma that hepatocytes increase the secretion of CRP (C-reactive protein) in response to myelomaderived cytokines, and the binding of CRP to FCGR2 on the surface of myeloma cells activates myeloma cells to promote osteoclastogenesis and bone destruction in vivo. It has also been shown that aggregation of FcgammaR on human bone marrow cells leads to tyrosine phosphorylation of inositol polyphosphate containing SHIP-2 (SRC homology 2 domain), which reveals a new regulatory role of the expression and function of inositol phosphatase SHIP-2 on FcgammaR mediated activation of human bone marrow cells [51]. These studies suggested that the FCGR2A we screened out might be highly associated with GIONFH and thus needs further study and validation. By consulting the literature, we could not find that the BNIP3L and ASXL1 were directly related to the occurrence of osteonecrosis. However, the protein encoded by BNIP3L/ NIX (NIX named Bcl-2/E1B19 kDa-interacting protein 3-like, which is based on 56% homology with BNIP3) is itself a pro-apoptotic protein. BNIP3L/NIX can not be detected in normal tissues, and they were found in a variety of organelle cultures in hypoxia-induced experiments [52,53]. The Bcl-2 gene family of hypoxiainduced NIX is expressed during erythropoiesis, and the researchers found that NIX is highly regulated during erythroid end-stage maturation [54]. This shows that NIX associated with BNIP3L can cause erythropenia under hypoxic conditions, which coincides with the ischemia and hypoxia that cause osteonecrosis of the femoral head. ASXL1 belongs to an oncogene with unstable genome and recent studies have shown that mutations in ASXL1 are found in hematopoietic cells of various myeloid tumors such as myelodysplastic syndrome and chronic myelomonocytic leukemia [55,56]. Our bold guess here is: are patients with mutations in ASXL1 more likely to develop GIONFH? While more evidence is needed. Therefore, we selected these four core genes with the highest correlation to GIONFH as research targets.

Fig. 4 According to
Through the induction of hormone and endotoxin, rabbit femoral head necrosis was caused. In the experiment later, we fully proved the success of animal model establishment through the diagnostic criteria  of magnetic resonance and pathological examination.
In the experimental group, we took four New Zealand white rabbits at different time points of the second, fourth and eighth weeks after hormone use for MRI examination, and dissect their femoral head specimens for pathological examination and gene detection.
Through the experiment, we were pleased to find that the quantitative expression of the four core genes, compared with the blank control group, in the animal samples of steroid-induced necrosis of the femoral head, the expression of ASXL1 and BNP3L is reduced, while the expression of FCGR2A and TYROBP is increased and hypolipidemic drugs can play a certain role in prevention and treatment of steroid-induced necrosis of femoral head in model animals. This result not only confirms that the core genes screened in the earlier stage do have an important relationship with the research purpose, and it is very important for early diagnosis with steroid-induced necrosis of femoral head. The prevalence of GIONFH is increasing, and people are diagnosed with it at a much earlier age. With the loss of function of the hip joint, it has brought great